clear all;
t = 0:1:100;
numSamples = 100;
nu = 0.3583e-4;
lambda = -nu*1^2;
sigma = sqrt(2.1983e-2);
Alpha = zeros(length(t),numSamples);
Alpha2 = zeros(length(t),numSamples);
Alpha(1,:) = 0;
Alpha2(1,:) = 0;
for j = 1:numSamples
    fprintf(1,'#%d sample\n',j);
    for i = 2:length(t)
        Alpha(i,j) = solveAlpha(Alpha(i-1,j),lambda,sigma,[t(i-1),t(i)],0.01);
        Alpha2(i,j) = solveAlpha2(Alpha2(i-1,j),lambda,sigma,[t(i-1),t(i)],0.01);
    end
end
Alpha2_a = Alpha.^2;

save solveAlpha_data
%%
figure(1);
mean_Alpha2_a = mean(Alpha2_a,2);
mean_Alpha2_b = mean(Alpha2,2);
plot(t,mean_Alpha2_a,'LineWidth',2);
set(gca,'FontSize',20);
xlabel('time(s)');
ylabel('mean(\alpha_1^2)')
% legend('From \alpha model','From \alpha^2 model');
print('fig_mean_Alpha2','-dpng');

%%
figure(2);
subplot(2,1,1);
hist(Alpha2_a(end,:));
xlabel('mean(\alpha_1^2), from \alpha model');
subplot(2,1,2);
hist(Alpha2(end,:));
xlabel('mean(\alpha_1^2), from \alpha^2 model');
print('hist_Alpha2','-dpng');

%%
figure(3);
var_Alpha2_a = var(Alpha2_a,ones(1,numSamples),2);
var_Alpha2_b = var(Alpha2,ones(1,numSamples),2);
plot(t,var_Alpha2_a,'LineWidth',2);
set(gca,'FontSize',20);
xlabel('time(s)');
ylabel('var(\alpha_1^2)')
% legend('From \alpha model','From \alpha^2 model');
print('fig_var_Alpha2','-dpng');